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Inspired by the collective phenomenon of territorial emergence, whereby animals move and inter- 
act through the scent marks they deposit, we study the dynamics of a ID Brownian walker in a 
random environment consisting of confining boundaries that are themselves diffusing anomalously. 
We show how to reduce, in certain parameter regimes, the non-Markovian, many-body problem of 
territoriality to the analytically tractable one-body problem studied here. The mean square dis- 
placement (MSD) of the ID Brownian walker within subdiffusing boundaries is calculated exactly 
and generalizes well known results when the boundaries are immobile. Furthermore, under certain 
conditions, if the boundary dynamics are strongly subdiffusive, we show the appearance of an in- 
' teresting non-monotonicity in the time dependence of the MSD, giving rise to transient negative 

, diffusion. 

^ . PACS numbers: 05.10.Gg, 05.40.Fb, 89.20.-a 

I. INTRODUCTION AND MOTIVATION 

' How disorder affects the intrinsic diffusion of a particle is a question of extreme practical relevance. It is often 
Oh, encountered in transport processes [ij whenever one attempts to confront theoretical predictions with actual experi- 
(-H ■ ments. Not surprisingly, the study of disorder has attracted a lot of interest in a variety of contexts and often with the 
aim of understanding the validity of the theoretical assumptions of idealized systems 0-111 • Although many different 
kinds exist, disordered systems are broadly divided into two classes depending on the type of heterogeneity: static 
or dynamic Q. When the heterogeneity is due to the presence of some 'quenched' obscrvablcs that do not change 
over time, one talks about a system with static disorder (see e.g. Q)- On the other hand, when the heterogeneity 
' emerges as the system evolves in time, one talks about dynamic disorder (see e.g. Q)- Of interest here is one of 
^ , the simplest dynamically disordered system: a random walker constrained to move within an area whose extent is 
' randomly fluctuating. 

I A problem where dynamic disorder plays a fundamental role is the formation of territorial patterns in the animal 
I kingdom. A recent study 0] shows how territories (spatial inhomogeneity) collectively emerge from the spatio- 
temporal trajectories of individual animals with interactions between them mediated via scent marks. In that study 
a microscopic stochastic model of territorial random walkers is introduced whereby each animal performs a random 
■ walk with the added ingredient that, when an animal visits a site, it leaves its scent that lasts for a fixed time 
Tas- During this time, no other animal can visit this site. As a result, each animal has a territory at any given instant 
t. This territory is the set of sites visited by the animal in the time period [t — T/jSjt]. An animal can diffuse only over 
sites that belong to its own territory and those sites which do not belong to the territory of any other animal. Thus, 
both the animals and their respective territories move with time, but their growth rates are different. The ccntroids 
of the territories subdiffusc, due to the exclusion principle |lll4l3{ . while the animals themselves move diffusively. 
. One of our aims is to reduce, in the ID case, the many-body non-Markovian problem of territory formation 
to the one-body problem of a random walker moving within a random environment. The wildly different time 
scales over which territories and animals displace, allows us to find the appropriate one-body effective description. 
If we focus on a specific individual, we can view it as being subject to reflecting boundary conditions, located 
at the position of the left and right neighbouring territorial borders, that are changing gradually relative to the 
time scale over which the walker diffuses. This gradual change in the external conditions (environment) is what 
characterizes an adiabatic process. We thus invoke an adiabatic approximation, which consists of taking the probability 
distribution of the walker to be identical to the one governed by a diffusion equation within an interval equal to the 
instantaneous separation distance between the neighboring territorial borders. The joint probability distribution 
P{x, Li, L2,t) of the walker position x and, respectively, the left and right boundary locations Li and L2 can be 
written as P{x, Li, L2,t) « Q{Li, L2,t)W{x,t\Li, L2), with W{x,t\Li, L2) being the walker probability distribution 
for a given value of Li and L2, and with Q{Li, L2,t), the joint distribution of the two boundary positions, to be 
determined. Our goal here is to show that it is possible to capture the dynamics of the random walkers in the 
territoriality problem by representing the dynamics of (5(Li,L2jO through a Fokker-Planck formalism [l4| with 
time-dependent diffusion coefficients. 

The paper is organized as follows. The analytic expression for the probability P{x, Li, L2, t) of the simplified one- 
body problem is computed in Sec. |lTl The results in Sec. |lT]are compared, in Sec. IIIIl with stochastic simulations of 
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two territorial random walkers in an interval with periodic boundary conditions. In Sec. llVl the analytic computation 
of the MSD is developed to show the appearance of transient negative diffusivity under certain conditions. The paper 
ends with a brief discussion in Sec. |Vl 



II. JOINT PROBABILITY DISTRIBUTION OF WALKER AND BOUNDARIES 



Since territorial boundaries undergo an exclusion process, the probability distribution of the separation distance 
between the left and right boundaries can be represented at long times by a Gaussian whose width is increasing 
proportionally to \Jt rather than or by a diffusion equation with an appropriate time-dependent diffusion constant 
[iSl [l6j . The mean of this Gaussian is centered around the average territory size, which is simply the inverse of 
the population density [§]. In our reduced one-body problem we mimic the presence of a mean territory size and 
fluctuations around this mean by assuming that the left and right boundaries are attached by a spring, whose rest 
position equals the mean territory size L. This leads us to the use of a Fokker-Planck formalism with a time-dependent 
diffusion constant and with a constraining quadratic potential UiL^^L-i) ~ ^{t){Li — Li — i)^/4 of the form 

The 2D force that J7(Li, L2) generates, that is ^{t){L2 — ii — L)/2 along Li and —j{t){L2 — Li — L)/2 along L2, has 
the effect of making the mean distance between the two boundaries equal L at long times. Counter to this driving 
force is the boundary spreading due to the random fluctuations associated with the positive function (p{t). For the 
time dependence of j{t) we focus on the case in which 7(f) = "ff(t) and to the constant case, i.e. when "f{t) = 7. In 
the former case we talk about a subordinated stochastic process since Eq. ([T]) can be rewritten via the transformation 
T = ds (p{s) as the evolution over the time t of a Fokker-Planck equation without time-dependent coefficients. In 
the latter case we have a time scale disparity between the diffusive process and the time it takes for the boundaries 
to reach their asymptotic separation value L. This, we will show in Sec. IIVI has the effect of generating negative 
diffusivity under certain conditions. 

Through the variable transformation X = L2 — L1 and £ = (i2 + ii)/2, Eq. ^ can be decoupled, i.e. Q{Li,L2, t) = 
Qi (£, t) Q2 (A, t) , into a differential equation for the boundary separation A and a differential equation for the boundary 
centroid location C. These two equations, together with the condition that, at all times, the two boundaries cannot 
cross each other, i.e. VQ{Li, L2,t) • n = where n is the normal to the line of points where Li = L2 can be 
solved exactly (see details in Appendix A). For localized initial conditions of Dirac 5 type of the form Q{Li, L2, 0) = 
S{Li + L/2)S{L2 — L/2) with £2,0 = L/2 — — Li.o, the time-dependent solution of Eq. ([T]) with the mentioned 
boundary condition is given by 

e b(t^ i_ g b(t) g c(t) 
Q[\C.t)=E[X) _, (2) 

where h(t) = %K J* ds(/3(s)e-2(G(t)-G(s)) ^-^^^i G{t) = /J ds7(s), c{t) = 2K Lp{s)ds and the Heaviside function H{y) 
is such that H{y) = I {H{y) = 0) if 1/ > (if y < 0). From the form of Eq. ^ it is evident that b{t) controls the 
diffusion of the boundary separation around the mean value L, and c{t) regulates the diffusion of the boundary centroid. 

For the choice ^{t) = ^ip{t) we have b{t) = A{K/^) |l — cxp —27 /J ds (p{s) |, whereas b{t) = 8K ds e~^^(*~''V(s) 
when 7(t) — 7. 

Armed with the probability distribution Q{Li, L2,t) one can now write down the full probability P{x, Li, L2,t) 
because W{x, t\Li, L2) can be obtained analytically with the method of images [l^ [2^ as 

W,,{x,t) = [H{x-L^)-H{x~L2)] J2 — +e ^ 

where xq is the initial walker position and w{t) = ADt with D being the diffusion coefficient of the walker. The 
Heaviside functions in Eq. ^ show explicitly that outside the region Li < x < L2 the walker probability distribution 
is identically zero. For ease of notation we have omitted here and in the rest of the paper Li, L2 in the definition of 
W . The quantity directly comparable to the simulation output, which is also easily accessible from field data, is the 
marginal probability distribution M[x,t) of the walker position x 
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M(x,t)= \ d\Z{\t) \ dL ^=W^,(x,t), (4) 
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wherein 

Z(A, t) = -^=2 cosh \-rh\- (5) 

Wc study in particular the dynamics of Eq. Q for a time dependence in the diffusion constant of the form 
(p{t) — a{t/()"~^ with C being a characteristic time. This choice aUows the tuning of the anomalous diffusion of 
the boundaries in terms of only one parameter, the exponent a. From the expression derived after performing the 
integration over C (see Appendix B for details) one notices that by using the dimensionless parameters D' = D/{L'^^), 
K' = K/{L'^j), 13 = and r = jt, the time-dependence in Eq. Q is lumped into the functions c{t)/L?, b{t)/L'^ and 
w{t)/L'^. D' and K' represent the average area covered in a time by, respectively, the walker and the boundaries 
relative to the square of the average boundary separation L. (3, being proportional to 7, represents the dimensionless 
rate at which the boundary separation returns to its average value. For subdiffusive processes, which we focus on here, 
we have < a < 1, with a — 1 being the limiting diffusive case. We thus have b{t)/L^ — 4A'' [l — exp(— 2/3^~"t")] 
when 7(t) = jip{t) or b{t)/L'^ = 8K' jS^'" dp exp[-2(T - p)]ap"-i when -f{t) = 7, whereas c{t)/L'^ = 2K'f3^-°'T°' 
and for the walker we have w{t)/L'^ = AD't. 




FIG. 1: Dependence of the parameters in the reduced model of Sec. |TT]in terms of the output from the stochastic simulations of 
the many-body problem. Panel (a) is a log-linear plot of the dimensionless K/D versus the dimensionless quantity Z = T^gp'^ 
where p' is the population density multiplied by the lattice spacing and T^g is the dimensionless quantity obtained by multiplying 
Tas with the random walk transfer rate F between nearest neighbour sites. Panel (b) is a log-linear plot of the dimensionless 
quantity s* /L^ versus Z^^'^ . The solid lines are least square fits of all the data points for respectively hog{K/D) = 1.11 — 0.95 
in panel (a) and Log(s/L'^) = 0.64 — 2.32 x Z^^'^ in panel (b). Each point in (a) is found by averaging over 10''' stochastic 
realisations of the simulation for sufficiently many time-steps to reach the boundary MSB's asymptotic regime. The time to 
reach this asymptotic regime increases exponentially with both p' and T'^g. For lower values of p', a high T'^g is required for 
the boundary MSD to reach the asymptotic regime before the limits of the box-size are reached. Therefore using either very 
high or very low p' is not practically possible, allowing us only the limited set of values in this plot. For (b), a single simulation 
of 10* time-steps gave sufficient data for finding each data-point. 



III. STOCHASTIC SIMULATIONS OF THE MANY-BODY PROBLEM 



We compare stochastic simulations of two territorial random walkers on a ring lattice with the reduced model of 
Sec. |TT]in the subordinate case, i.e. when ^[t) — j(p{t) = ^■yt^^ , With this choice of j{t), the MSD of the boundary 
separation in the simulations is qualitatively similar to the reduced model, in that ((L — A)^) is monotonic, increasing 
initially and saturating at long times. In order to make this comparison we need to determine the dependence of 
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K and 7 on the active scent time Tas and the walker population density p, which in this case is simply twice the 
inverse of the number of lattice sites in the ring. In Fig. [T] we show these dependencies for a range of values of 
active scent time and population density by plotting the dimensionless parameter K/D and the asymptotic value 
s* = s{t — > +oo)/L^ of the normalized MSD of the boundary separation distance s{t) (see Appendix A for its analytic 
expression). These parameters depend on the quantity Z = T^gp'^ where T^g and p' are dimensionless versions of 
Tas and p respectively (see caption of Fig. [1]). Qualitatively one would expect K to decrease as function of Z as 
the latter represents the time for the scent mark to disappear, Tas- divided by the mean first passage time to move 
between the left and right boundaries, which is proportional to a first approximation to the length squared of the 
domain itself [2l| . The rate of movement of the boundaries is obviously larger the longer it takes the walker to refresh 
its own scent marks as well as the shorter the scent remains active. The fitting lines obtained in Fig. [T] provide 
the means for selecting the appropriate parameter values K' , D' and /3, which are present in the expression for the 
marginal probability distribution of the walker in Eq. (jl9p . 
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FIG. 2: Comparison of the walker probability distribution from the simulated model (dashed lines) with the reduced model 
(solid lines) at time tF=10,500 for four values of Z: a) Z=1.8, b) Z=1.4, c) Z=1.2 and d) Z=0.6 To plot the probability 
distribution of the reduced model, M{x,t) in Eq. [19] is multiplied by the lattice constant. In all panels the density has been 
kept fixed at p' = 0.02, whereas T^g = 9000,7000,6000,3000 for the panel a), b), c) and d), respectively. Running the 
simulations requires a choice of initial scent spatial profile, i.e. the time value associated with each lattice site that defines how 
long before the scent stops being active. We have selected a biologically relevant scent profile with minimum values at the scent 
boundaries and a maximum in the middle of the territory. The initial location of the two walkers is consequently the middle 
of their respective territories (see ref. for further details). 

For a choice of T^^ and p' in the simulations, one determines K/{-^L'^) with L = p^^, and K/D of the reduced 
model. The walker diffusion constant is simply o?F where a is the lattice spacing and F the random walk transfer 
rate between nearest neighbour sites. F and a do not need to be measured since they are input of the model. From 
Fig. [1^ one can now determine K and subsequently 7 from Fig. [T]d. In Eq. ([!]) for convenience K has units of a 
diffusion constant and ip{t) is a dimensionless quantity since time is rescaled with the time constant C,. To measure K 
from the MSD of the boundary in the simulation, it is necessary to choose a value for C,. Figure [1^ was produced by 
picking C — 1 / F so this value of C must be used in any comparisons between the simulation output and the reduced 
model. However, this choice is arbitrary since if we were to choose C, — C /F for some C > instead, each value of K 
measured from the simulation output would be reduced by a factor of C^~", i.e. points in Fig. [T^ would be shifted 
down by Log(C^~"£'). Then, if the new values of C and K were placed back into Eq. ([1]), the changes would cancel 
one another out. 

In Fig. [2] we show the result of this exercise for four sets of T'^g — p' parameter choice by superimposing the 
simulated probability distribution of either walker and the marginal probability distribution of the reduced model, 
whose parameters have been selected from the empirically derived long-time relations between K/D and Z and s* jl? 
and Z. By looking sequentially at the panels from (a) to (d) in Fig. [5] it is evident that the dissimilarity between the 
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theoretical curves and their respective simulated curves increases as Z gets smaller. This is due to the assumptions 
in our reduced model. The adiabatic approximation relies upon the walker probability distribution to be close to the 
case when boundaries are immobile, a situation corresponding to Tas = oo. A reduction of the value of Z corresponds 
to a decrease in the active scent time and/or a decrease of the walker population density. For a given density p of 
walkers, a large reduction in Tas causes the region where foreign scent is present to move too quickly compared to the 
diffusion time of the walker. Scent deposited at a boundary site is more likely to become inactive before the animal 
has had a chance to re-visit that site, causing the boundary to move with higher probability. Similarly, for a given 
value of Tasi ^ big reduction in p makes the walkers move within a much larger area. The time it takes for the walker 
to roam between the territorial boundaries becomes too large to assume that Li and L2 move little compared to the 
walker. In other words, as Z diminishes the conditions for the assumptions in the adiabatic approximation breaks 
down. 



IV. NON-MONOTONIC DEPENDENCE OF THE MEAN SQUARE DISPLACEMENT 

An important quantity that characterizes the walker anomalous movement is the time-dependence of the MSD 
averaged over all boundaries' locations: 

dL2 / dL, / dx{x^xofQ{Li,L2,t)W,,ix,t), (6) 

-00 J — 00 J Li 

where the symbol (((...))) indicates that the average is over the three variables x, Li and L2. By performing two of 
the integrations Eq. (jS]) can be reduced (see Appendix C for details) to the analytic formula 

((((.-cm - ^l + J^ + ii + ^ + l ^AZ(A,t)g|-^cos(^^je ^ 

,r''(2n-l)''[c(t) + .^(t)l 1 

which is a generalization of the time-dependence of a Brownian walker MSD within immobile boundaries whose 
distance equals L. In the limit if — >■ we have in fact that c{t) — > and b{t) — > which transforms Z{X,t) into 
(5(A — i), reducing Eq. (O to a well known expression (see e.g. ref. j22|). 

An inspection of Eq. ([7]) when xq = shows that by expanding the exponential up to second order in the even 
series and up to first order in the odd series one obtains (((x^))) ^ w{t) for t 0. The short-time dependence of the 
MSD is thus dependent only on the walker movement when away from the boundaries. The long time behaviour on 
the other hand shows that (((a;^))) ~ + 6(t)/24 -I- c{t)/2. The two series in Eq. ([7]) in fact decay to zero when 

t +00 because of the exponential terms. In other words the boundary dispersal dictates how the MSD increases at 
long times through the diffusion of the boundary separation b{t) and the boundary centroid c(t). Since b{t) and c{t) 
contains (p{t), the choice of the exponent a allows us to control the time-dependence of the (normalized) asymptotic 
factor a{t) = 1 + b{t) / (2i^) + 6c{t) / and consequently the long-time anomalous diffusion characteristics of the MSD. 
For the case 7(<) = 7<p(i), b{t) reaches a positive asymptotic value, as shown previously, whereas b{t) approaches at 
long times when "f{t) = 7. This means that in both cases the dominant effect in the walker MSD at long times is due 
to the boundary centroid rather than the boundary separation. In the remaining part of this section we focus on the 
case 7(t) =7. 

The non-monotonicity of the time-dependence in the b{t) under certain conditions may also appear in a{t). When 
that occurs we say that the walker undergoes the phenomenon of negative diffusion as a result of the interplay between 
two contributions to the boundary movement: the centroid displacement, governed by c{t), and the deviation of the 
boundary separation distance from L, governed by b{t). The former is monotonically increasing, whereas the latter 
is not, increasing initially, before reaching a peak and decaying to zero at long times. If at any point the increase of 
c{t) is sufficiently slow compared to the decrease of b{t), the boundary MSD will decrease. 

The parameter dependence for this non-monotonicity can be found by first studying the time dependence of the 
walker's MSD asymptotic expression a(<), which possesses a maximum and a minimum whenever the a-dependent 
function Vq{t) > 2 with Uq(t) = r^~" JJ" ds exp[— 2(t — s)]s"~^. A numerical study of Va{T) shows that this occurs 
when a < ao with ao ~ 0.105 and irrespective of any other parameters. In Fig. [3]we show a few examples with a < ao 
where the MSD possesses a dip in time for different values of the parameters K' , D' , /? and a. Fig. [3^ shows that, 
for fixed 13', /? and a, decreasing K' tends to smooth out the appearance of the dip, eventually making it disappear 
with further decrease of K' (not shown in the figure). On the other hand, maintaining K' , /? and a fixed, it is the 
increase of D' that eventually makes the dip disappear, as one can see by looking sequentially at the various curve in 
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FIG. 3: Non-monotonicity of the walker's MSD (Eq. ((7]) with xo — 0) for different parameter values. The subfigure (a) shows 
how the non-monotonic behaviour of the MSD, for two choices of a < ao, tends to be washed away as K' decreases, while the 
D' and /? are kept fixed. In the subfigure (b) the MSD dependence over D' is shown while keeping the other parameters fixed: 
K' = 0.01, Q = 0.05, 13 = 0.1, and in subfigure (c) the MSD dependence over /3 is shown with K' = 0.01, a = 0.05, D' = 1. 



Fig. ^jp. A qualitative similar dependence is depicted in Fig. [5J; as /? decreases while keeping the other parameters 
fixed. 

To explain why varying if', D' and /? as shown in Fig. [3] make the transient negative difFusivity present in the 
expression a(t) disappear, we rewrite Eq. ([7]). We do so by performing the A integration (see Appendix D for details). 
In the resulting expression besides the asymptotic term a{t), the time dependence is lumped into powers of order 
f{t)^/'^ with fc > 1, but more importantly terms like exp[— L^/6(t)], e^xp[—2^T^Jf{t)] and exp[— tt^/ f{t)] containing 
the dimensionless quantity f{t) = [c{t) + w{t)\/h{t) = (r" + 2D'tI3"-^ / K') / {4 /J" dp exp[-2(T - p)]ap°'-'^}. Since 
f{t) is linearly proportional to D' and inversely proportional to K' (3^~", whereas L'^/b{t) is inversely proportional to 
K'j3^~°', it is possible to understand qualitatively the direct and inverse proportionality on K' and D' that can be 
evinced by looking sequentially at the various curves in Fig. |31 As D' increases, the value of f{t) increases, and the 
MSD more rapidly becomes very close to its asymptotic expression a{t), as clearly shown by the trend in the plots 
depicted in Fig. [3^. The opposite trend is observed when i^T'/?^^" increases due to its inverse proportionality in the 
MSD expression. 

Having identified the parameter combination for which the MSD time-dependence is non-monotonic, we now try 
to interpret it, and in particular we try to explain the reasons for the sudden transition at a < ao. For that, it 
is sufficient to consider the asymptotic expression a{t)L^ /12 of the MSD and determine from it the value of the 
(time-dependent) diffusion constant. With simple algebra one can rewrite d/ dt[a{t) / 12] as an effective diffusion 
coefficient Deff{t) = i^(/>(t){l + T[a]Ei^ai-2'Yt)/3} where E^^siz) = ^21=0 l^l^^n b) is the 2-parameter Mittag- 
Leffler function. K(f>(t), the first term in Deff{t), is the time-dependent diffusion coefficient of either boundary, whereas 
the second term in Deff{t) is the one for the boundary separation. Interestingly, the function T[a]Ei^a{—'2jt)/3 starts 
at 1/3 at t = and remains positive for any values of time only if a = 1, i.e. when the Mittag-Leffler reduces to an 
exponential. For a < 1 on the other hand, T[a]Ei,a{—2'-ft) /3 becomes negative, reaches a minimum and eventually 
approaches the value zero from below as r(a)[6r(a — 1)7^]"^. 

The negative diffusion is due to the qualitative difference in time scales of two processes: the random movement 
of the boundary separation variable A and its rate of return to the value L. The former is sublinear, except when 
a = 1, whereas the latter is always linear. At short times the Gaussian nature of the process makes the effective 
diffusion constant of the boundary separation positive. Although at time i = we have set the random variable A to 
be (5-localized, i.e. zero everywhere except at A = L, at any time t > 0, A can assume values much larger than L due 
to the Gaussian shape of Q2{X,t). Irrespective of the shape of the constraining potential, the diffusion constant for 
the boundary separation is initially positive making the MSD increase sub-linearly (when a < 1). The spring force 
acts against this initial transient increase in the MSD by pushing the value A back to L. The drift towards L occurs 
on such a fast time scale that it counteracts the boundary separation dispersion, making the MSD decrease its value 
after a certain amount of time. At that moment the system can be described by an effective diffusion constant with 
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FIG. 4: Lines in parameter space D' — a, for different values of K'j3^ delineating the regions below which the walker's MSD 
is non-monotonic. 



a negative sign. 

Although the above description is valid for any a < 1, the transition is observed only at a = ao rather than a = 1. 
This is because we have only described the dynamics of the boundary separation, neglecting the effects of the boundary 
centroid. The decreasing trend observed in the full asymptotic expression a{t) eventually ends once the reduction in 
the MSD of the boundary separation is overcome by the monotonic increase in the MSD of the boundary centroid. 
The effective diffusion constant Deff{t) is in fact the sum of an always positive contribution due to the boundary 
centroid; and a contribution due to the boundary separation, which is initially positive and negative afterwards. When 
one considers both contributions, it turns out that only when the boundary subdiffusion is particularly pronounced, 
i.e. when a < ao ~ 0.105, does the negative diffusion due to the spring force fail to compensate for the positive 
diffusion of the boundary centroid, making the MSD non-monotonic. 



V. CONCLUSIONS 



We have analyzed the movement of a ID Brownian walker within anomalously diffusing boundaries, focusing on 
the case in which each boundary moves subdiffusively at all times. We have shown how the analytic dependence of 
the walker MSD depends on the movement statistics of the boundary separation and boundary centroid by using a 
Fokker-Planck formalism with a time-dependent diffusion coefficient. The study has been motivated by the need to 
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derive an approximate model of the many-body phenomenon of territorial emergence consisting of walkers interacting 
through a renormalized exclusion process for which individuals are excluded from an area recently visited by other 
walkers Q. 

Despite the non-Markovian nature of territorial dynamics, we have exploited the time scale disparity between the 
territory movement, which is subdiffusive, and the walker movement, which is diffusive, to construct a reduced one- 
body model that captures the salient features of the many-body problem. Through an adiabatic approximation we 
have shown how the analytic model of the single walker moving within subdiffusing boundaries is able to represent the 
complex dynamics present in animal territoriality. A comparison between the probability distribution of the analytic 
model and the stochastic simulations of two territorial random walkers on a ring has confirmed the validity of the 
adiabatic approximation. 

Explicit expressions of the walker MSD have been computed, which are direct generalizations of the MSD of random 
walker within fixed boundaries. An interesting outcome of our formalism is the appearance, under certain conditions, 
of a non-monotonic behaviour of the Brownian walker MSD at intermediate times, that is to say, the system undergoes 
the phenomenon of transient negative diffusion (seen in other context e.g. [13 )■ 

Although a direct application of our equations is the displacement of an animal within its own moving territorial 
boundaries any situation where the movement of two or more objects are separated by an enclosing interface or 
boundaries (see e.g. could benefit from our study. We have focused on the case of a Brownian walker within 

subdiffusing boundaries, but our formalism could also be employed in other applications in which the boundary 
movement is modelled as a subordinate anomalous diffusion process, by choosing more general dependence of the 
boundary time-dependent diffusion constant ip{t). Similarly, one could take into account anomalously diffusing walkers 
by having w{t) = AD dsxis) where Dx{t) is a time-dependent diffusion constant. 

Applications of these kinds may occur at much smaller spatial scales compared to those typical of animal territorial 
dynamics. At the sub-cellular level for example, the positioning of certain cell components is crucial to its functioning 
and relies upon the continuous exploration through a variety of motor proteins of the slowly changing intracellular space 
[2^. The moving boundaries in this context would represent the translation and deformation of the cell membrane. 
Examples at the cellular level can be found in developmental studies on certain macrophages that move via contact- 
inhibition interaction, a form of repulsion mechanism upon contact. This mechanism allows the macrophages to 
partition the space they continuously monitor for the presence of infectious agents [1^. This space partitioning 
creates effective boundaries between macrophages that continuously change in time, which could also be modeled 
though the help of the formalism developed here. 

We have chosen to represent anomalous diffusion via a subordinate Gaussian process, but other representations are 
possible, one such being a time non-local formalism given by the Generalized Master Equation (GME) or its 

equivalent, the Fractional Diffusion Equation (FDE) [2^ (see e.g. ref. [s^l for the equivalence between GME and 
FDE). A comparative investigation between the time local formalism used here, i.e. dQ{z,t)/dt = K(p{t)W'^Q{z,t), 
and the time non-local GME formalism, i.e. dQ{z,t)/dt = K ds(j>{t — s)V^Q(z, s), where (/)(t) is the so-called 
memory function, has been carried out in ref. [3l|. That comparative analysis, however, did not consider the presence 
of a constraining potential. We plan to extend that work by performing a comparative analysis between a time 
local and a time non-local representation of the boundary fluctuations in presence of a constraining force, as well as 
analyzing the influence that the potential shape has on the results presented here. 
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Appendix A 

The variable transformation to the boundary separation A and boundary ccntroid C described in Sec. HIl allows one 
to rewrite Eq. ^ as 



and the corresponding initial conditions as Q(£,A, 0) = d{C — CQ)d{X — Aq) with Cq = (^i.o + i2.o)/2 and Aq = 
L2.0 — il.O- 
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Eq. (|8]) can now be decoupled as Q{Li, L2,t) = Qi{C,t)Q2{X,t), into a differential equation for the boundary 
separation A 



dQiiCt) _ K ...d^QiiC^ 
dt " 2 ' 

and a differential equation for the boundary centroid location C 

Eq. © and (|10p arc solved together with the reflection condition at A = 0, i.e. 

d 



dX 



Q2(A,t)|A=0 -0. 



(9) 



(10) 



(11) 



The solution of Eq. ([9]) trivially gives a Gaussian whose width increases in time proportionally to c(t), whereas the 
solution of Eq. (|10p is more involved as it needs to be solved with the method of characteristics [S^l- The general 
solution of Eq. (fTO|) is given by 



Q2(A,t) 



(A-A(t))^ 

e Ht) 



(12) 



where 



X{t) = L + e-^^'HXo-L), 



h{t) 



8K 



dsv?(s)e-2(G(*)-G(.)) 



(13) 
(14) 



and G{t) is defined so that G'{t) = 7(t). Once again with the help of the method of images [l^l one takes into account 
the reflective boundary condition at A = and obtains the general solution of Eq. ([T0|) and (fTTj) as 



(A-A(t))^ 



Q2iX,t)^H{X)- 



+ e 



(A + A(f))^ 



(15) 



Eq. ([T5|) represents the probability distribution of the positive distance A and it is given by the sum of two Gaussian 
curves with the same width but centered around two different values: L for the first one and —L for the second one 
(for the choice of initial conditions we have that X{t) = L). If ^(t) = ■y(p{t) the width of each of these Gaussian 
curves starts at and then increases monotonically to the value AK/^ if ^{t) — jip{t). On the other hand when 
j{t) ~ 7, these curve widths have an initial increase to a maximum followed by a decrease back to 0. Given that the 
second Gaussian is centered at some negative value of L, for values of L that are not too small one can think of the 
dynamics of Q2(A,t) as given essentially by the spreading of the Gaussian centered around the positive constant L. 
The deviations from this qualitative description can be determined exactly by computing the MSD s{t) = ((A — i)^) 
of A from L. A simple integration gives 



sit) 




"<*) +4-— - 

b{t) Kt) 



1 - erf I 



L 



(16) 



which shows that the MSD equals b{t)/2, as one would expect for the case of a Gaussian with variance b{t)/2 and 
mean L, only when L/ \/b{t) >> 1. When comparing with simulations, in the case "/{t) = jip{t), the pertinent value 
is the limit s(t — > +oo) = s* . To find this limit, simply substitute AK/'y for b{t) everywhere in Eq. ([T6)) . 

Appendix B 

The marginal probability distribution of the walker can be computed via the double integral expression 



M{x,t) = / dXZ{X,t) 
Jo 



+ 00 



dC 



e <=<*) 



g,, {x, t) [H{x -C + A/2) -H{x-C- A/2)] 



(17) 
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where gxo{x,t) represents the series terms in Eq. ([3]). By making use of the Poisson summation formula [l^l for the 
first series in gx„{x,t) one can transform Eq. (jl7p into 



M(x,t) = / d\Z{\t) 







.-A/2 2A 



e 4^-^ 



+ 00 

E 



n— — oo 



[x^2£+(2,i+l)A + xol'' 



(18) 



Integration over L then gives 



dAZ(A,t) 







erf 



+ 00 



a; — xo 
A 



+ 00 

E 



[i + iO + (2n+l)A]^ 
g 4c(t) + m(t) 



erf 



h{x,xo,t) + Xq„ (t) 



with 



erf 



h{x,xo,t)-X (^q+{t)+4^) 



^Ac{t) + w{t) 



(19) 



Hx,xo,t)=xlM + 2M^ 



2xn 



Icitl 
wit)'' 



1 cit) 



2 y c(i) y w{t) 



(20) 



Appendix C 



For a simphfied analytic expression of Eq. ([6]) , it is convenient to rewrite each infinite series in Eq. ([3]) by making 
use of the Poisson summation formula [l^l , which gives 



1 1 ^ 
a + aE 



71=1 



cos nn 



X - Xq 



A 



+ (-l)"cos mr 



X + Xo 



c 



cos 2mr- 



+(-l)"sin nn 



.To 



A 



sm I 2n7r— 



(21) 



To obtain Eq. (O from ^ in the main text, one can perform the integration over x between the values Li and 
L2, and subsequently the £- integration after transforming to (A,£)-space. Alternatively one can exploit the relation 

((x-xo)2) = \^-S^T{w,„{x,t)} + 2ixo-i^T{Wx,{x,t)} + xlT{W,,{x,t)}\ \k=Q, where T{Wx,{x,t)} is the Fourier 

transform of Wxo{x,t). We write here the details of this calculation due to its potential applicability when only the 
Fourier expression of Wxq {x, t) is known, e.g. when the walker probability is a Levy distribution. 

Since W^^ [x, t) is the product of g^^ (t, t) and the step function H{x ~ Li) — H{x — L2), the Fourier expression of 
Wxf,{x,t) is the convolution oi T{Wo{x,t)} with F{H{x — Li) — H{x — L2)} which are equal, respectively, to 



mr 

T 



I e" ~ 5\k 

n=l 



6 k- 



e >- S \ k + 



)]}' 



and 



T{H{x - Li) - H{x - L2)} (k) = 



2e''=^ sin(^^ 



By exploiting the presence of the various Dirac delta function in 
domain between Eq. and and obtain 



T{Wx,{x,t)}{k) 



2e*''^sin(^) 1 
kX ^ A 



E 



^I^, (22) 
(23) 

one can compute the convolution in Fourier 
,--^A'~^)^sin[{k-^-f)^] 
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AY 



e'(''+^)^sin[(fc + ^) A]' 
fc+ if 



From the Fourier expression F{Wxo{x,t)} in it is now straightforward to determine that 

2n7r(£ — .tq) 



A2 +2? 4A2(-1)" 



+ 00 



^ (2n)27r2 

n— 1 ^ ^ 



COS 



A 



-8A(/: - xo) V — -2-^ sin 



(2n- l)7r(£-a;o) 
A 



(2,1-1)-' ^{t) 



(24) 



(25) 



After inserting Eq. (|25p in Eq. ([6|), changing the Li- and i2-integration into a A- and £- integration and performing 
the £-integration, one then recovers Eq. ^ in the main text. 



Appendix D 

When xo = the MSD expression ([7]) can be simphfied further if the infinite summation and the integral can be 
interchanged, which occurs when the series is uniformly convergent [33| (this occurs for all times when ^{t) = ^ip{t), 
whereas it occurs for t > when ^{t) = 7). One first needs to notice that 



+00 



dX A^"e 



=2 (V™ /e^^"V'3 



(26) 



When Eq. (pS)) is evaluated at /3 = 1 one obtains 

f + OO 



dXe-^ -I^A^" = ^-n-e™(2a)e-^", 



(27) 



where 



2™ ^ (to - kV.kl 



(28) 



are Bessel polynomials [3J]. By using the series representation of the hyperbolic cosine one can now perform the 
A-integration and obtain 



7n—0 



7r2 



hit) ±2? ©rn+l [2nn^f{t)j (_i)„ _ 



g-27rnV/(t)_|_ 



4c(t)^ Q-i-(2»-l)x// ft)j (-1)" 



(2n-l)V7(*) 



2r7. - 1 



(29) 



with fit) = [c(t) + w(t)]/6(t) as defined in the main text. The summation over n can now be performed by rewriting 
the two series as, respectively, 



+ 00 +00 /'_g-2ir^/(t)A 

^"T7'=-2„-27rny7{t) _ V / 



^(-l)"n^-2e-2W/(*) = ^ 



n=l 



n=l 



n 



2-fe 



(30) 



and 



^(-l)"(2n - i)fe-ie— (2"-i)%/iW = ^2''-\--^/m 



i-fc 



(31) 
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and noticing that '^^=1 z"'/n^ = Lir{z) is the polylogarithm [3^, which reduces to a polynomial for any integer r < 0, 
and that X^n^^ -^"/ + '^Y ~ ^ {z,r,a) is the so-called Lerch trascendent function [s^ , which reduces to a ratio of y 
polynomials for any r < when a = 1/2. From Eq. ([501 and ([5l]) it follows that 



^ 12+^ + ^ + 4!;^^ (^^fen 



2-fc)! 



TT ""^VKOy (2m)! ^ (m-fc)!fc! 



fc)!fc! 



$ -e 



1 



1 - /fc, ■ 



) 

(32) 



To simplify further the infinite series expression in p2p one can sum the infinite m-series for specific values of k. For 
example, if all the fc = and fc = 1 terms in the first series and all the fc = and fc = 1 terms in the second series are 
summed over all m- values, Eq. ([5^ becomes 



1 

12 



Li, -e-^-Vm 



In ( 1 + er'^Vm 



b{t) 



1 

24 



Li, f-e-W/W 



27r2 



b{t) ^ 



c(0 

m+l 



1 4 

arctan 

2 TT 



47r2 



/_L^Y"' 1 (2m + 2-fc)! 



E 



\\h{t)) (2m)! ^ (m + 1 - fc)!fc! 



4^V7(I)] Li2-, (-e-2-x/7w) 



2c(t) f L 



(2m- fc)! 

-'\b(t)) (2?7i)! ^ (m- fc)!fc! 



E 



E 



4^y7W e-"V/(*)$ -e-2-\//(t),i-fc 



(33) 



Continuing this process for all fc gives the following single-sum expression, involving the generalized hypergeometric 
function 2-^2(01, 02; 61, &2; z) (see e.g. [s^l and references therein or [s!] for the properties of the 2-F2 hypergeomteric 
function) . 



1 

12 



Li, (-e-2V/(*) 



ln(l + e-2'rv/(*: 



1 

24 



Li2 -e-2V/(*) 



27r2 



pin(l + e-2-x/7w) 



-c{t) 



1 4 

arctan 

2 TT 

2\ 



(e-"V7(*)^ - 4V7(t) (^1 - e"^^ - 



+ g-27rV/(t) 



fe=2 



(2fc-2)!V b{t) 



'E 

fc=2 



1 + -, i + -; fc, fc - i; 4^ki2-,.(-e-2-V7W) 
^^^ 2'2 2'' 2' 6(t)/ ^ ' 

2\ 



(2fc)! V b{t) 



F2U + il + r^k + l,k+l;^)H-e-'-^\l~k^-) (34) 



2' 2 2 



2' b{t) 



[1] M. Dresden, Rev. Mod. Phys. 33, 265 (1961). 

[2] S. Alexander, J. Bernasconi, W. R. Schneider, and R. Orbach, Rev. Mod. Phys. 53, 175 (1981). 
[3] J. W. Haus and K. W. Kehr, Phys. Rep. 150, 236 (1987). 

[4] H. E. Stanley and N. Ostrowsky (eds.) Random Fluctuations and Pattern Growth: Experiments and Models (Kluwer 
Academic, Dordreclit, 1988). 

[5] V. May and O. Kiihn, Charge and Energy Transfer Dynamics in Molecular Systems (Wiley- VCH, Berlin, 2000). 



13 



[6] V. M. Kenkre, in Exciton Dynamics in Molecular Crystals and Aggregates: the Master Equation Approach, edited by G. 

Hohler, Springer Tracts in Modern Pliysics, Vol. 94, (Springer, Berlin, 1982), pp. . 
[7] B. H. Hughes, Random Walks and Random Environments, Vol. 2 (Oxford Univ. Press, 1996). 

[8] A. R. Bishop, D. K. Campbell and S. Pnevmatikos, Disorder and Nonlinearity Springer Proceedings in Physics Vol. 39 
(Springer, Berlin, 1989). 

[9] L. Giuggioh, J. R. Potts, S. Harris, PLoS Comp. Biol. 7(3): el002008. doi: 10. 1371/journal.pcbi. 1002008 (2011). 
[10] A. Okubo and S. A. Levin, Diffusion and Ecological Problems: Modern Perspectives 2nd edition, (Springer, 2002). 
[11] T. E. Harris, J. Appl. Probab. 2, 323 (1965). 
[12] C. Landim, Ann. Probab. 20, 206 (1992). 

[13] G. Schonherr and G.M. Schiitz, J. Phys. A: Math. Gen. 37, 8215 (2004). 

[14] H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer, Berlin, 1989). 

[15] G. K. Batchelor, Proc. Cambr. Philos. Soc. 48 (1952) 345. 

[16] B. B. Mandelbrot and J. W. Van Ness, SIAM Rev. 10, 422 (1968). 

[17] L. Giuggioh, G. M. Viswanathan, V. M. Kenkre, R. R. Parmenter and T. L. Yates, Europhys. Lett. 77, 40004 (2007). 
[18] T. Ambjornsson, L. Lizana, M. A. Lomholt and R. J. Silbey, J. Chem. Phys. 129, 185106 (2008). 
[19] S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943). 

[20] E. W. MontroU, B. J. West, in Fluctuation phenomena edited by E.W. MontroU and J.L. Lebowitz (Amsterdam, North- 
Holland, 1987), pp. 61-206. 
[21] S. Condamin, O. Benichou, M. Moreau, Phys. Rev. E 75, 021111 (2007). 

[22] L. Giuggioh, G. Abramson, V. M. Kenkre, G. Suzan, E. Marce and T. L. Yates, BuU. Math. Biol. 67, 1135 (2005). 
[23] V. G. Karpov, Phys. Rev. Lett. 75, 2702 (1995). 

[24] D. R. Nelson, T. Piran and S. Weinberg (eds.). Statistical Mechanics of Membranes and Interfaces, 2nd ed., (World 

Scientific, Singapore, 1989). 
[25] 1. M. Tolic-N0rrelykke, Curr. Opin. CeU. Biol. 22, 21 (2010). 

[26] B. Stramer, S. Moreira, T. Millard, L Evans, C.-Y. Huang, O. Sabet, M. Milner, G. Dunn, P. Martin and W. Wood, J. 
CeU. Biol. 189, 681 (2010). 

[27] V. M. Kenkre in Statistical Mechanics and Statistical Methods in Theory and Application, edited by U. Landman (New 

York, Plenum, September 1977), pp. 441-461. 
[28] V. M. Kenkre, in Modem Challenges in Statistical Mechanics: Patterns, Noise, and the Interplay of Nonlinearity and 

Complexity, edited by V. M. Kenkre and K. Lindenberg (New York, AlP, 2003), pp. 63-100. 
[29] R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000). 

[30] L. Giuggioh, F. Sevilla and V. M. Kenkre, J. Phys. A: Math. Theor. 42, 434004 (2009). 

[31] V. M. Kenkre and F. J. Sevilla, in Contributions to Mathematical Physics: a Tribute to Gerard G. Emch, edited by T. S. 

Ah and K. B. Sinha, (Hindustani Book Agency, New Delhi 2007) pp. 147-160. 
[32] P. H. Moon and D.E. Spencer, Partial Differential Equations (Heath, Lexington, 1969). 
[33] K. Knopp, Theory and Applications of Infinite Series (New York, Dover, 1990). 
[34] E. Grosswald, Bessel Polynomials (Lecture Notes in Mathematics) (New York, Springer 1978). 
[35] A. Jonquiere, B. Soc. Math. Fr. 17, 142 (1889). 
[36] M. Lerch, Acta Math. 11, 19 (1887). 
[37] R. B. Paris, J. Comput. Appl. Math. 173, 379 (2005). 

[38] Hardy, G. H. in Ramanujan: Twelve Lectures on Subjects Suggested by His Life and Work, 3rd ed. (New York, Chelsea 
1999) pp. 101-112. 



